%% config
rootdir = ['..' filesep '211007_Faruq_2FLIBats'];

right2left = 1;

fudgeFactor = 0.95;

manual = 0;
manualPos = [202 1080];
%% find sample directories
if rootdir(end) ~= filesep
    rootdir(end+1) = filesep;
end
if exist([rootdir 'Online' filesep],'dir')
    fovdir{1} = [rootdir 'Online' filesep];
else
    listDir = dir(rootdir);
    fovdir = cell(length(listDir)-2,1);
    for k = 3:length(listDir)
        if exist([rootdir listDir(k).name filesep 'Online' filesep],'dir')
            fovdir{k-2} = [rootdir listDir(k).name filesep 'Online' filesep];
        end
    end
    fovdir = fovdir(~cellfun('isempty',fovdir));
end

%% load
for idxDir = 1:length(fovdir)
fprintf('%s\n',fovdir{idxDir});
listFoVBegin = dir([fovdir{idxDir} 'M*_FoV_begin.png']);
listFoVEnd = dir([fovdir{idxDir} 'M*_FoV_end.png']);
listFileNames = [{listFoVBegin.name} {listFoVEnd.name}]';
measNum = nan(size(listFileNames));
for k = 1:length(listFileNames)
    strName = listFileNames{k};
    idx1 = 1;
    idx2 = find(strName == '_');
    measNum(k) = str2double(strName(idx1+1:idx2(1)-1));
end
measNum = unique(measNum);
for idxMeas = 1:length(measNum)
fprintf('M%d\n',measNum(idxMeas));
fprintf('type\tdate\t\t\tposition inlet\tposition outlet\tchannel legth\n');
file1 = dir([fovdir{idxDir} sprintf('M%d_FoV_begin.png',measNum(idxMeas))]);
if ~isempty(file1)
    img1 = double(imread([file1.folder filesep file1.name]));
    try
        chPos1 = funcFindChannelPosition(1, img1, 1, [], fudgeFactor);
    catch
        chPos1 = [1 size(img1,2)];
    end
    title('Begin');
    fprintf(['begin\t' file1.date '\t' sprintf('%d',chPos1(2)) ' pix' '\t' sprintf('%d',chPos1(1)) ' pix' '\t\t' sprintf('%d',diff(chPos1)) ' pix = ' sprintf('%g',diff(chPos1)*0.34) ' �m\n']);
end
fileMeas = dir([fovdir{idxDir} sprintf('M%d_*.tdms',measNum(idxMeas))]);
if ~isempty(fileMeas)
    fprintf(['meas ' '\t' fileMeas.date '\n']);
end
file2 = dir([fovdir{idxDir} sprintf('M%d_FoV_end.png',measNum(idxMeas))]);
if ~isempty(file2)
    img2 = double(imread([file2.folder filesep file2.name]));
    try
        chPos2 = funcFindChannelPosition(2, img2, 1, [], fudgeFactor);
    catch
        chPos2 = [1 size(img2,2)];
    end
    title('End');
    fprintf(['end  \t' file2.date '\t' sprintf('%d',chPos2(2)) ' pix' '\t' sprintf('%d',chPos2(1)) ' pix' '\t\t' sprintf('%d',diff(chPos2)) ' pix = ' sprintf('%g',diff(chPos2)*0.34) ' �m\n']);
end

if ~isempty(file1) && ~isempty(file2)
    repeatSelection = 1;
    while repeatSelection
        fprintf('change:\n');
        fprintf('x1: %d pix -> %d pix = %g �m\n',chPos1(1),chPos2(1),(chPos2(1) - chPos1(1)) * 0.34);
        fprintf('x2: %d pix -> %d pix = %g �m\n',chPos1(2),chPos2(2),(chPos2(2) - chPos1(2)) * 0.34);
        fprintf('[0] find positions manually\n[1] begin\n[2] end\n[3/enter] mean\n');
        choice = input('select dataset: ');
        if choice == 0
            choiceManual = input('manual mode: determine for [1] begin, [2] end: ');
            if any(choiceManual == 1)
                manualPos = input('begin, position [x1, x2]: ');
                chPos1 = funcFindChannelPosition(1, img1, 1, manualPos, fudgeFactor);
                title('Begin');
            end
            if any(choiceManual == 2)
                manualPos = input('end, position [x1, x2]: ');
                chPos2 = funcFindChannelPosition(2, img2, 1, manualPos, fudgeFactor);
                title('End');
            end

            fprintf(['begin\t' file1.date '\t' sprintf('%d',chPos1(2)) ' pix' '\t' sprintf('%d',chPos1(1)) ' pix' '\t\t' sprintf('%d',diff(chPos1)) ' pix = ' sprintf('%g',diff(chPos1)*0.34) ' �m\n']);
            fprintf(['end  \t' file2.date '\t' sprintf('%d',chPos2(2)) ' pix' '\t' sprintf('%d',chPos2(1)) ' pix' '\t\t' sprintf('%d',diff(chPos2)) ' pix = ' sprintf('%g',diff(chPos2)*0.34) ' �m\n']);

        elseif choice == 1
            chPos = chPos1;
            repeatSelection = 0;
        elseif choice == 2
            chPos = chPos2;
            repeatSelection = 0;
        else
            chPos = round(0.5*(chPos1 + chPos2));
            repeatSelection = 0;
        end
    end
    save([fovdir{idxDir} sprintf('M%d_',measNum(idxMeas)) 'channelPos.mat'], 'chPos');
end
if ~isempty(file1) && isempty(file2)
    chPos = chPos1;
    save([fovdir{idxDir} sprintf('M%d_',measNum(idxMeas)) 'channelPos.mat'], 'chPos');
end
if isempty(file1) && ~isempty(file2)
    chPos = chPos2;
    save([fovdir{idxDir} sprintf('M%d_',measNum(idxMeas)) 'channelPos.mat'], 'chPos');
end

end
end

return
%% functions
function chPos = funcFindChannelPosition(fig, img, right2left, manualPos, fudgeFactor)
%% load image
figure(fig);clf;
imagesc(img);axis image;
colormap gray;
%% gradient
if isempty(manualPos)
[gradX, gradY] = gradient(img);
imgGrad = sqrt(gradX.^2 + gradY.^2);
%% edge detection
[~,threshold] = edge(imgGrad,'sobel');
BWs = edge(imgGrad,'sobel',threshold * fudgeFactor);
% imshow(BWs)
%% mask margins
lmargin = round(0.03*size(BWs,2));
rmargin = lmargin;
BWs(:,1:lmargin) = 0;
BWs(:,size(BWs,2)-(rmargin-1:-1:0)) = 0;
% imshow(BWs)
%% dilate the image
se90 = strel('line',3,90);
se0 = strel('line',3,0);
BWsdil = imdilate(BWs,[se90 se0]);
% imshow(BWsdil)
%% fill interior gaps
BWdfill = imfill(BWsdil,'holes');
% imshow(BWdfill)
%% remove small objects
BW1 = bwareaopen(BWdfill, 100);
BW1Area = regionprops(BW1,'Area');
BW1Area = sort([BW1Area.Area]);
if length(BW1Area) >= 3
    BW2 = bwareaopen(BWdfill, round(mean(BW1Area(end-2:end-1))));
else
    BW2 = BW1;
end
% imshow(BW2)
%% 
imgMod = BW2;
sizeImg = size(img);
%% upper half
upperHalf = imgMod(1:round(sizeImg(1)/2),:);
% imshow(upperHalf);
upperY = sum(upperHalf .* repmat((1:size(upperHalf,1))',1,size(upperHalf,2)),1) ./ sum(upperHalf,1);
% plot(upperY);

upperYGrad = gradient( meanfilt1( upperY, 2*round(0.04*length(upperY) / 2)+1 ) );
riseIdxs = find(upperYGrad > 0.4);
fallIdxs = find(upperYGrad < -0.4);
% hold on;
% yLim = get(gca,'YLim');
% line(riseIdxs(1)*[1 1],yLim,'Color','r');
% line(riseIdxs(end)*[1 1],yLim,'Color','r');
% line(fallIdxs(1)*[1 1],yLim,'Color','r');
% line(fallIdxs(end)*[1 1],yLim,'Color','r');
% hold off;

try
    x = 1:sizeImg(2);
    fit1 = polyfit(x(riseIdxs(1):riseIdxs(end)),upperY(riseIdxs(1):riseIdxs(end)),1);
    fit2 = polyfit(x(riseIdxs(end):fallIdxs(1)),upperY(riseIdxs(end):fallIdxs(1)),1);
    fit3 = polyfit(x(fallIdxs(1):fallIdxs(end)),upperY(fallIdxs(1):fallIdxs(end)),1);

    x11 = -(fit1(2) - fit2(2)) / (fit1(1) - fit2(1));
    x12 = -(fit3(2) - fit2(2)) / (fit3(1) - fit2(1));

    set(0, 'currentfigure', fig); yLim = get(gca,'YLim'); xLim = get(gca,'XLim'); hold on;
    plot(x,fit1(1)*x+fit1(2),'r');
    plot(x,fit2(1)*x+fit2(2),'r');
    plot(x,fit3(1)*x+fit3(2),'r');
    hold off;
    set(gca,'YLim', yLim); set(gca,'XLim', xLim);
catch
    x11 = [];
    x12 = [];
end

%% lower half
yOffs = round(sizeImg(1)/2);
lowerHalf = imgMod(yOffs:end,:);
% imshow(lowerHalf);
lowerY = sum(lowerHalf .* repmat((1:size(lowerHalf,1))',1,size(lowerHalf,2)),1) ./ sum(lowerHalf,1);
% plot(lowerY);

lowerYGrad = gradient( meanfilt1( lowerY, 2*round(0.04*length(lowerY) / 2)+1 ) );
riseIdxs = find(lowerYGrad > 0.4);
fallIdxs = find(lowerYGrad < -0.4);
% hold on;
% yLim = get(gca,'YLim');
% line(riseIdxs(1)*[1 1],yLim,'Color','r');
% line(riseIdxs(end)*[1 1],yLim,'Color','r');
% line(fallIdxs(1)*[1 1],yLim,'Color','r');
% line(fallIdxs(end)*[1 1],yLim,'Color','r');
% hold off;

try
    x = 1:sizeImg(2);
    fit1 = polyfit(x(riseIdxs(1):riseIdxs(end)),lowerY(riseIdxs(1):riseIdxs(end)),1);
    fit2 = polyfit(x(fallIdxs(end):riseIdxs(1)),lowerY(fallIdxs(end):riseIdxs(1)),1);
    fit3 = polyfit(x(fallIdxs(1):fallIdxs(end)),lowerY(fallIdxs(1):fallIdxs(end)),1);

    x22 = -(fit1(2) - fit2(2)) / (fit1(1) - fit2(1));
    x21 = -(fit3(2) - fit2(2)) / (fit3(1) - fit2(1));

    set(0, 'currentfigure', fig); yLim = get(gca,'YLim'); xLim = get(gca,'XLim'); hold on;
    plot(x,fit1(1)*x+fit1(2) + yOffs-1,'r');
    plot(x,fit2(1)*x+fit2(2) + yOffs-1,'r');
    plot(x,fit3(1)*x+fit3(2) + yOffs-1,'r');
    hold off;
    set(gca,'YLim', yLim); set(gca,'XLim', xLim);
catch
    x22 = [];
    x21 = [];
end

end

%% mean inlet and outlet position
if isempty(manualPos)
    if isempty(x11) && not(isempty(x21))
        pos_inlet = x21;
    elseif not(isempty(x11)) && isempty(x21)
        pos_inlet = x11;
    else
        pos_inlet = 0.5*(x11 + x21);
    end
    
    if isempty(x12) && not(isempty(x22))
        pos_outlet = x22;
    elseif not(isempty(x12)) && isempty(x22)
        pos_outlet = x12;
    else
        pos_outlet = 0.5*(x12 + x22);
    end
else
    pos_inlet = min(manualPos);
    pos_outlet = max(manualPos);
end

set(0, 'currentfigure', fig); yLim = get(gca,'YLim');
hold on;
line(pos_inlet*[1 1],yLim,'Color','b');
line(pos_outlet*[1 1],yLim,'Color','b');
hold off;

%% result
if ~isempty(manualPos)
    pos_inlet = min(manualPos);
    pos_outlet = max(manualPos);
end

if right2left
    pos_inlet = size(img,2) + 1 - pos_inlet;
    pos_outlet = size(img,2) + 1 - pos_outlet;
end

% fprintf('position inlet = %g pix\n', pos_inlet);
% fprintf('position outlet = %g pix\n', pos_outlet);
% fprintf('channel length = %g pix = %g um\n', abs(pos_outlet-pos_inlet), abs(pos_outlet-pos_inlet)*0.34);

chPos = [round(min(pos_inlet,pos_outlet)) round(max(pos_inlet,pos_outlet))];
% save([rootdir filesep 'Online' filesep 'channelPos.mat'], 'channelPos');
end